function power_iteration()
% POWER_ITERATION 幂迭代法演示
% 
% 演示幂迭代法及其变形用于计算特征值和特征向量

fprintf('=== 幂迭代法 ===\n\n');

%% 1. 基本幂迭代法
fprintf('1. 基本幂迭代法\n');
fprintf('迭代公式: x_{k+1} = A*x_k / ||A*x_k||\n');
fprintf('收敛到最大模特征值对应的特征向量\n\n');

% 测试矩阵
A = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('矩阵 A =\n');
disp(A);

% 精确特征值（用于比较）
eigenvals_exact = sort(eig(A), 'descend');
fprintf('精确特征值: [%.6f, %.6f, %.6f]\n', eigenvals_exact);
fprintf('主特征值: λ₁ = %.6f\n', eigenvals_exact(1));

% 幂迭代
x0 = [1; 1; 1];  % 初始向量
max_iter = 20;
tol = 1e-10;

fprintf('\n幂迭代过程:\n');
fprintf('迭代    特征值近似    误差\n');
fprintf('----    ----------    ----\n');

x = x0 / norm(x0);
lambda_old = 0;

for k = 1:max_iter
    % 幂迭代步骤
    y = A * x;
    lambda = x' * y;  % Rayleigh商
    x = y / norm(y);  % 归一化
    
    % 计算误差
    error = abs(lambda - eigenvals_exact(1));
    
    if k <= 10 || mod(k, 5) == 0 || error < tol
        fprintf('%4d    %10.6f    %.2e\n', k, lambda, error);
    end
    
    % 收敛检查
    if abs(lambda - lambda_old) < tol
        fprintf('收敛于迭代 %d\n', k);
        break;
    end
    
    lambda_old = lambda;
end

fprintf('最终特征向量: [%.6f; %.6f; %.6f]\n', x);

%% 2. 收敛速度分析
fprintf('\n2. 收敛速度分析\n');
fprintf('收敛速度取决于 |λ₂/λ₁| 的比值\n\n');

ratio = abs(eigenvals_exact(2) / eigenvals_exact(1));
fprintf('|λ₂/λ₁| = %.6f\n', ratio);
fprintf('理论收敛率: %.6f\n', ratio);

% 分析实际收敛率
errors = [];
x = x0 / norm(x0);
for k = 1:15
    y = A * x;
    lambda = x' * y;
    x = y / norm(y);
    error = abs(lambda - eigenvals_exact(1));
    errors = [errors; error];
end

% 计算实际收敛率
if length(errors) > 2
    actual_rates = errors(2:end) ./ errors(1:end-1);
    avg_rate = mean(actual_rates(5:end));  % 忽略前几次迭代
    fprintf('实际平均收敛率: %.6f\n', avg_rate);
end

%% 3. 反幂迭代法
fprintf('\n3. 反幂迭代法\n');
fprintf('计算最小模特征值: (A - σI)^(-1) * x\n\n');

sigma = 0;  % 位移参数
A_shifted = A - sigma * eye(size(A));

fprintf('位移 σ = %.1f\n', sigma);
fprintf('反幂迭代计算最小特征值:\n');

x_inv = x0 / norm(x0);
fprintf('迭代    特征值近似    误差\n');
fprintf('----    ----------    ----\n');

for k = 1:10
    % 反幂迭代步骤
    y = A_shifted \ x_inv;  % 求解线性系统
    mu = x_inv' * y;        % Rayleigh商
    x_inv = y / norm(y);    % 归一化
    
    lambda_inv = 1/mu + sigma;  % 转换回原特征值
    error_inv = abs(lambda_inv - eigenvals_exact(3));
    
    fprintf('%4d    %10.6f    %.2e\n', k, lambda_inv, error_inv);
end

%% 4. 带位移的反幂迭代
fprintf('\n4. 带位移的反幂迭代\n');
fprintf('计算接近给定值的特征值\n\n');

target = 3.2;  % 目标特征值
sigma_shift = target;
A_shift = A - sigma_shift * eye(size(A));

fprintf('目标特征值: %.1f\n', target);
fprintf('位移 σ = %.1f\n', sigma_shift);

% 找到最接近目标的精确特征值
[~, closest_idx] = min(abs(eigenvals_exact - target));
closest_eigenval = eigenvals_exact(closest_idx);
fprintf('最接近的精确特征值: %.6f\n', closest_eigenval);

x_shift = x0 / norm(x0);
fprintf('迭代    特征值近似    误差\n');
fprintf('----    ----------    ----\n');

for k = 1:8
    y = A_shift \ x_shift;
    mu = x_shift' * y;
    x_shift = y / norm(y);
    
    lambda_shift = 1/mu + sigma_shift;
    error_shift = abs(lambda_shift - closest_eigenval);
    
    fprintf('%4d    %10.6f    %.2e\n', k, lambda_shift, error_shift);
end

%% 5. Rayleigh商迭代
fprintf('\n5. Rayleigh商迭代\n');
fprintf('自适应位移的反幂迭代，三次收敛\n\n');

x_ray = x0 / norm(x0);
lambda_ray = x_ray' * A * x_ray;  % 初始Rayleigh商

fprintf('Rayleigh商迭代:\n');
fprintf('迭代    特征值近似    误差\n');
fprintf('----    ----------    ----\n');

for k = 1:6
    % 使用当前特征值估计作为位移
    A_ray = A - lambda_ray * eye(size(A));
    
    try
        y = A_ray \ x_ray;
        x_ray = y / norm(y);
        lambda_ray = x_ray' * A * x_ray;  % 更新Rayleigh商
        
        % 找到最接近的精确特征值
        [~, closest_idx] = min(abs(eigenvals_exact - lambda_ray));
        error_ray = abs(lambda_ray - eigenvals_exact(closest_idx));
        
        fprintf('%4d    %10.6f    %.2e\n', k, lambda_ray, error_ray);
        
    catch
        fprintf('迭代 %d: 矩阵接近奇异\n', k);
        break;
    end
end

%% 6. 同时迭代法（子空间迭代）
fprintf('\n6. 同时迭代法\n');
fprintf('同时计算多个主特征值和特征向量\n\n');

% 计算前两个特征值
p = 2;  % 子空间维数
X0 = randn(size(A, 1), p);  % 初始子空间
[X, ~] = qr(X0, 0);  % 正交化

fprintf('同时计算前%d个特征值:\n', p);
fprintf('迭代    λ₁近似      λ₂近似      误差₁       误差₂\n');
fprintf('----    -------     -------     -----       -----\n');

for k = 1:12
    % 子空间迭代
    Y = A * X;
    [X, R] = qr(Y, 0);
    
    % 计算Rayleigh-Ritz近似
    H = X' * A * X;  % 投影矩阵
    [V_sub, D_sub] = eig(H);
    ritz_values = sort(diag(D_sub), 'descend');
    
    error1 = abs(ritz_values(1) - eigenvals_exact(1));
    error2 = abs(ritz_values(2) - eigenvals_exact(2));
    
    if k <= 8 || mod(k, 2) == 0
        fprintf('%4d    %7.4f     %7.4f     %.2e    %.2e\n', ...
                k, ritz_values(1), ritz_values(2), error1, error2);
    end
end

%% 7. 幂迭代的变形：Arnoldi迭代预览
fprintf('\n7. Arnoldi迭代预览\n');
fprintf('构造Krylov子空间的正交基\n\n');

% Arnoldi过程（简化版本）
m = 4;  % Krylov子空间维数
v1 = x0 / norm(x0);
V = zeros(size(A, 1), m);
H = zeros(m+1, m);
V(:, 1) = v1;

fprintf('Arnoldi过程构造 %d 维Krylov子空间:\n', m);

for j = 1:m
    w = A * V(:, j);
    
    % Gram-Schmidt正交化
    for i = 1:j
        H(i, j) = V(:, i)' * w;
        w = w - H(i, j) * V(:, i);
    end
    
    H(j+1, j) = norm(w);
    if j < m
        V(:, j+1) = w / H(j+1, j);
    end
    
    fprintf('步骤 %d: ||w|| = %.6f\n', j, H(j+1, j));
end

% 计算Ritz值
H_trunc = H(1:m, 1:m);
ritz_arnoldi = sort(eig(H_trunc), 'descend');
fprintf('Arnoldi Ritz值: [%.4f, %.4f, %.4f, %.4f]\n', ritz_arnoldi);
fprintf('前两个的误差: [%.2e, %.2e]\n', ...
        abs(ritz_arnoldi(1:2) - eigenvals_exact(1:2)));

%% 8. 收敛性分析和比较
fprintf('\n8. 收敛性分析和比较\n');

% 不同方法的收敛历史
methods = {'Power', 'Inverse', 'Rayleigh'};
convergence_data = cell(3, 1);

% 重新运行各方法记录收敛历史
% 幂迭代
x = x0 / norm(x0);
power_errors = [];
for k = 1:15
    y = A * x;
    lambda = x' * y;
    x = y / norm(y);
    power_errors = [power_errors; abs(lambda - eigenvals_exact(1))];
end
convergence_data{1} = power_errors;

% 反幂迭代
x = x0 / norm(x0);
inverse_errors = [];
for k = 1:15
    y = A \ x;
    mu = x' * y;
    x = y / norm(y);
    lambda = 1/mu;
    inverse_errors = [inverse_errors; abs(lambda - eigenvals_exact(3))];
end
convergence_data{2} = inverse_errors;

% Rayleigh商迭代（重新运行）
x = x0 / norm(x0);
lambda = x' * A * x;
rayleigh_errors = [];
for k = 1:8
    try
        A_temp = A - lambda * eye(size(A));
        y = A_temp \ x;
        x = y / norm(y);
        lambda = x' * A * x;
        [~, idx] = min(abs(eigenvals_exact - lambda));
        rayleigh_errors = [rayleigh_errors; abs(lambda - eigenvals_exact(idx))];
    catch
        break;
    end
end
convergence_data{3} = rayleigh_errors;

% 绘制收敛曲线
figure('Name', '幂迭代方法收敛比较');
colors = {'b-o', 'r-s', 'g-^'};
for i = 1:3
    if ~isempty(convergence_data{i})
        semilogy(1:length(convergence_data{i}), convergence_data{i}, ...
                 colors{i}, 'LineWidth', 2, 'MarkerSize', 6);
        hold on;
    end
end
xlabel('迭代次数');
ylabel('误差 (对数尺度)');
title('幂迭代方法收敛比较');
legend(methods, 'Location', 'best');
grid on;

%% 9. 应用：PageRank算法
fprintf('\n9. 应用：PageRank算法\n');
fprintf('Google PageRank算法的简化版本\n\n');

% 简单网络图的邻接矩阵
% 节点1 -> 节点2,3
% 节点2 -> 节点3
% 节点3 -> 节点1
adj_matrix = [0, 1, 1; 0, 0, 1; 1, 0, 0];
n_pages = size(adj_matrix, 1);

% 构造Google矩阵
damping = 0.85;
out_degrees = sum(adj_matrix, 2);
out_degrees(out_degrees == 0) = 1;  % 避免除零

G = damping * (adj_matrix ./ out_degrees) + (1-damping)/n_pages * ones(n_pages);
fprintf('Google矩阵 G =\n');
disp(G);

% 使用幂迭代计算PageRank
x_pr = ones(n_pages, 1) / n_pages;  % 初始均匀分布
fprintf('PageRank迭代:\n');
fprintf('迭代    页面1      页面2      页面3\n');
fprintf('----    -----      -----      -----\n');

for k = 1:10
    x_pr = G * x_pr;
    x_pr = x_pr / sum(x_pr);  % 归一化（保持概率分布）
    
    if k <= 6 || k == 10
        fprintf('%4d    %.4f     %.4f     %.4f\n', k, x_pr);
    end
end

fprintf('最终PageRank分数: [%.4f, %.4f, %.4f]\n', x_pr);

% 验证：这应该是主特征值为1的特征向量
eigenvals_G = eig(G);
fprintf('Google矩阵的主特征值: %.6f\n', max(real(eigenvals_G)));

end